统计计量 | 用R做多元线性回归分析（文末有福利）

回归分析是一种统计学上分析数据的方法，目的在于了解两个或多个变量间是否相关、相关方向与强度，并建立数学模型，以便通过观察特定变量（自变量），来预测研究者感兴趣的变量（因变量）。

多元线性回归的一般表现形式是：

其中，i为1,2，...，n；k为解释变量的数目，称为回归参数。

Q

回顾一下计量课上，多元线性回归模型的基本假定有哪些？

假设1：解释变量是非随机的或固定的，且各X之间互不相关（无多重共线性）；

假设2：随机误差项具有零均值、同方差及不序列相关；

假设3：解释变量与随机项不相关；

假设4：随机项满足正态分布；

假设5：样本容量趋于无穷时，各解释变量的方差趋于有界常数；

假设6：回归模型的设定是正确的。

如果违背以上假设，那么我们回归的统计结果很可能就不正确了。

那么利用R软件，如何进行正确而高效的多元线性回归分析呢？本期我们用R中的lm（）函数对个人抽样数据（有效样本1219个）进行拟合线性回归，样例数据展示如下表：

id	income	yrsch	hhwktot	age	gender	hanzu	urban	province
430001	1700	12	24.5	27	0	1	0	43
500002	3500	6	56	30	0	1	0	50
650003	4000	9	28	28	0	1	0	65
360004	5500	8	2.33	31	0	1	0	36
......	......	......	......	......	......	......	......	......
131215	1000	5	42	27	1	1	0	13
511216	2600	6	28	38	1	1	0	51
411217	1200	12	21	37	1	1	0	41
431218	6000	9	21	26	1	1	1	43
461219	4000	5	31.5	35	1	1	1	46

左右滑动查看更多

注：id是个人身份代码，income是个人收入，yrsch是教育年限，hhwktot是在家休闲时间，age是年龄，gender是性别，hanzu是民族虚拟变量，urban反映户籍制度，province是指省份。

基于样本数据，假如我们想探究个人收入和其他因素的关系，包括教育程度、休闲时间、年龄、性别、民族和户籍。现在我们使用lm（）函数拟合多元线性回归模型。

1.加载包并读取数据

## 加载包
library(openxlsx)
library(dplyr)
library(stargazer)
library(car)
## 读取数据
pid_data<-read.xlsx("OLS_Data.xlsx")

左右滑动查看更多

2.检验二变量关系

## 提取关键变量
corr_data<-as.data.frame(pid_data[,c("income","yrsch","hhwktot","age","gender","hanzu","urban")])
## 散点矩阵
scatterplotMatrix(corr_data,spread=F,lty.smooth=2)

左右滑动查看更多

图1 散点矩阵

scatterplotMatrix（）函数在非对角线区域绘制变量散点图，并添加平滑拟合线，对角线区域绘制每个变量的密度图和轴须图。从图1中可以看到，大多数连续变量是正态分布曲线，但都存在一定程度的偏斜。个人收入随着教育程度和年龄的增加而增加，随着闲暇时间的增加而减少等等。

3.OLS回归：

## lm（）函数多元回归
fit<-lm(log(income)~yrsch+hhwktot+age+gender+hanzu+urban,data=pid_data)
## 查看回归结果
summary(fit)

#Call:
#lm(formula = log(income) ~ yrsch + hhwktot + age + gender + hanzu +
#    urban, data = pid_data)

#Residuals:
#    Min    1Q    Median    3Q      Max
#-3.5591 -0.4155  0.0625  0.4997  2.0844

#Coefficients:
#             Estimate Std. Error  t value  Pr(>|t|)    
#(Intercept)  7.413813   0.209533  35.383   < 2e-16 ***
#yrsch        0.068765   0.008535   8.056  1.87e-15 ***
#hhwktot     -0.012416   0.001625  -7.640  4.39e-14 ***
#age          0.011887   0.005514   2.156  0.031297 *  
#gender      -0.172216   0.051564  -3.340  0.000864 ***
#hanzu       -0.067964   0.085215  -0.798  0.425283    
#urban        0.572335   0.051353  11.145   < 2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Residual standard error: 0.7517 on 1212 degrees of freedom
#Multiple R-squared:   0.35,Adjusted R-squared:  0.3468
#F-statistic: 108.8 on 6 and 1212 DF,  p-value: < 2.2e-16

## 导出回归结果，回归结果见图2
stargazer(fit,type="text",
          title="Regression Results",
          dep.var.labels=c("log(income)"),
          covariate.labels=c("yrsch","hhwktot",'age',
                             "gender","hanzu","urban"),
          ci=FALSE, ci.level=0.90, single.row=TRUE)

左右滑动查看更多

图2 回归结果

从图2看，教育年限的回归系数为0.069，且在1%水平上显著，表示控制闲暇时间、年龄、性别、民族、户口不变时，教育程度每增加一年，个人收入将会提升6.9%。此外，F统计量在1%水平上显著，调整R方为0.347，总体来看，所有预测变量解释了收入水平34.68%的方差。

4.模型诊断

上面的拟合回归模型是不是好的估计呢？我们知道，一个好的回归模型应当满足假设前提，包括：自变量要与因变量有线性关系；残差基本呈正态分布；残差方差基本不变（同方差性）；残差（样本）间相互独立。下面我们用plot（）绘图观察：

## 回归诊断
par(mfrow=c(2,2))
plot(fit)

左右滑动查看更多

图3 回归诊断

线性：“残差图与拟合图（Residuals vs Fitted）”反映了因变量和自变量线性相关性，从图3左上看，两者存在一个比较清晰的线性关系。

正态性：“正态Q-Q图（Normal Q-Q）”反映了当因变量成正态分布，则残差值也应该是一个均值为0的正态分布。正态Q-Q图（图3右上）是在正态分布对应值下，标准化残差的概率图，若满足正态分布，那么图中的点应该落在呈45度角的直线上。

同方差性：若满足不变方差的假设，那么在位置尺度图（Scale-Location，图3左下）中，围绕曲线周围的点应该随机分布。

独立性：我们无法从图3中分辨出因变量值是否相互独立，只能从收集的数据进行验证。

你可能会注意到，“残差与杠杆图（Residuals vs Leverage，图3右下）”会向我们传达哪些信息呢？其实，它反映出样本的离群点、高杠杆值点和强影响点。从图3右下观察，存在着这样的离群点。

相对plot（fit）方法，R中还有更好的R语言包去拟合和评价回归模型的能力，如cat、gvlma等。

_____________________________

文末福利：点击阅读原文，下载样本数据Excel。

►一周热文

数据Seminar

这里是大数据、分析技术与学术研究的三叉路口

作者：简华（何年华）审阅：杨奇明编辑：青酱

    欢迎扫描👇二维码添加关注